Examples

In [1]:

from bcn.binding_network import binding_network
from bcn.catalysis_network import catalysis_network, binding_and_catalysis

from bcn.binding_network import binding_network from bcn.catalysis_network import catalysis_network, binding_and_catalysis

In [2]:

import numpy as np
import pandas as pd
import time
import sympy as sp
from scipy.integrate import solve_ivp
from scipy.linalg import null_space
import itertools
import cvxpy as cp

import plotly.graph_objects as go
import plotly.express as px

import numpy as np import pandas as pd import time import sympy as sp from scipy.integrate import solve_ivp from scipy.linalg import null_space import itertools import cvxpy as cp import plotly.graph_objects as go import plotly.express as px

Define binding networks¶

There are three reaction networks need to idenfify

1. $E + S \leftrightharpoons C$
2. $E + 2S \leftrightharpoons C$
3. $E + S \leftrightharpoons C_1,\ S + C_1 \leftrightharpoons C_2$

In [3]:

def binding_catalysis_network(k_bind, k_cat, bn, total_E, total_const_idx, xcat_in_total_idx, cat_active_in_xbind_idx):
    """
    Using Binding and Catalysis Networks (BCN) to smiluate the speed of reaction
    inputs:
    outputs:
    """
    
    s_mat = np.array([
        [-1], # S
        [1], # P
    ])
    
    cn = catalysis_network(s_mat)

    kbind = np.array(k_bind) # kes (um)
    kcat = np.array(k_cat) # kcat_s2p
    
    bcn = binding_and_catalysis(bn, cn, kbind, kcat, np.array([total_E,]), total_const_idx, xcat_in_total_idx, cat_active_in_xbind_idx)
    
    return bcn

def binding_catalysis_network(k_bind, k_cat, bn, total_E, total_const_idx, xcat_in_total_idx, cat_active_in_xbind_idx): """ Using Binding and Catalysis Networks (BCN) to smiluate the speed of reaction inputs: outputs: """ s_mat = np.array([ [-1], # S [1], # P ]) cn = catalysis_network(s_mat) kbind = np.array(k_bind) # kes (um) kcat = np.array(k_cat) # kcat_s2p bcn = binding_and_catalysis(bn, cn, kbind, kcat, np.array([total_E,]), total_const_idx, xcat_in_total_idx, cat_active_in_xbind_idx) return bcn

define binding network: $E + S \leftrightharpoons C$¶

In [4]:

# binding part for 1st reaction network
# E + S = C

n_mat=np.array([
    [1, 1, 0, -1]  # binding reaction E + S <-> C
])

# x_sym = sp.symbols('x_E, x_S, x_P, x_C,')
# t_sym = sp.symbols('q_{E}, q_{S}, t_{P},')
# k_sym = sp.symbols('K_{ES},')

# bn1 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym)
bn1 = binding_network(n_mat, is_atomic=True)
bn1.l_mat

# binding part for 1st reaction network # E + S = C n_mat=np.array([ [1, 1, 0, -1] # binding reaction E + S <-> C ]) # x_sym = sp.symbols('x_E, x_S, x_P, x_C,') # t_sym = sp.symbols('q_{E}, q_{S}, t_{P},') # k_sym = sp.symbols('K_{ES},') # bn1 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym) bn1 = binding_network(n_mat, is_atomic=True) bn1.l_mat

Out[4]:

array([[ 1.,  0.,  0.,  1.],
       [ 0.,  1.,  0.,  1.],
       [ 0.,  0.,  1., -0.]])

define binding network: $E + 2S \leftrightharpoons C$¶

In [5]:

# binding part for 2nd reaction network
# E + 2S = C

# binding part

n_mat=np.array([
    [1, 2, 0, -1]  # binding reaction E + 2S <-> C
])

# x_sym = sp.symbols('x_E, x_S, x_P, x_C,')
# t_sym = sp.symbols('t_{E}, t_{S}, t_{P},')
# k_sym = sp.symbols('K_{ES},')
# bn2 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym)

bn2 = binding_network(n_mat, is_atomic=True)
bn2.l_mat

# binding part for 2nd reaction network # E + 2S = C # binding part n_mat=np.array([ [1, 2, 0, -1] # binding reaction E + 2S <-> C ]) # x_sym = sp.symbols('x_E, x_S, x_P, x_C,') # t_sym = sp.symbols('t_{E}, t_{S}, t_{P},') # k_sym = sp.symbols('K_{ES},') # bn2 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym) bn2 = binding_network(n_mat, is_atomic=True) bn2.l_mat

Out[5]:

array([[ 1.,  0.,  0.,  1.],
       [ 0.,  1.,  0.,  2.],
       [ 0.,  0.,  1., -0.]])

define binding network: $E + S \leftrightharpoons C_1,\quad S + C_1 \leftrightharpoons C_2$¶

In [6]:

# binding part for 3rd reaction network
# E + S = C1, S + C1 = C2

n_mat = np.array([
    [1, 1, 0, -1, 0], 
    [0, 1, 0, 1, -1],
])

# x_sym = sp.symbols('x_E, x_S, x_P, x_{C1}, x_{C2},')
# t_sym = sp.symbols('q_{E}, q_{S}, q_{P},')
# k_sym = sp.symbols('K_{ES}, K_{SC1},')
# bn3 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym)

bn3 = binding_network(n_mat, is_atomic=True)
bn3.l_mat

# binding part for 3rd reaction network # E + S = C1, S + C1 = C2 n_mat = np.array([ [1, 1, 0, -1, 0], [0, 1, 0, 1, -1], ]) # x_sym = sp.symbols('x_E, x_S, x_P, x_{C1}, x_{C2},') # t_sym = sp.symbols('q_{E}, q_{S}, q_{P},') # k_sym = sp.symbols('K_{ES}, K_{SC1},') # bn3 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym) bn3 = binding_network(n_mat, is_atomic=True) bn3.l_mat

Out[6]:

array([[ 1.,  0.,  0.,  1.,  1.],
       [ 0.,  1.,  0.,  1.,  2.],
       [ 0.,  0.,  1., -0., -0.]])

Define catalysis networks¶

define catalysis network: $S \rightarrow P$¶

In [7]:

s_mat = np.array([
    [-1], # S
    [1],  # P
])

cn = catalysis_network(s_mat)

s_mat = np.array([ [-1], # S [1], # P ]) cn = catalysis_network(s_mat)

Define binding and catalysis networks¶

In [8]:

total_E = 5  # total enzyme
total_S = 5 + 1e-6
total_P = 1e-6

total_E = 5 # total enzyme total_S = 5 + 1e-6 total_P = 1e-6

define BCN1¶

bn1 and cn

In [9]:

k_bind_1 = np.array([1])  # kbind = k- / k+
k_cat_1 = np.array([1])  # kcat_s2p
# bcn

total_const_idx_1 = np.array([0,])  # qe is first in totals (qe, qs, qp)
xcat_in_total_idx_1 = np.array([1, 2])  # qs is 2nd element in totals, tp is not contained.
cat_active_in_xbind_idx_1 = np.array([3])  # c_es

k_bind_1 = np.array([1]) # kbind = k- / k+ k_cat_1 = np.array([1]) # kcat_s2p # bcn total_const_idx_1 = np.array([0,]) # qe is first in totals (qe, qs, qp) xcat_in_total_idx_1 = np.array([1, 2]) # qs is 2nd element in totals, tp is not contained. cat_active_in_xbind_idx_1 = np.array([3]) # c_es

In [10]:

bcn1 = binding_and_catalysis(
    bn1,
    cn,
    k_bind_1,
    k_cat_1,
    np.array([total_E,]),
    total_const_idx_1,
    xcat_in_total_idx_1,
    cat_active_in_xbind_idx_1
)

bcn1 = binding_and_catalysis( bn1, cn, k_bind_1, k_cat_1, np.array([total_E,]), total_const_idx_1, xcat_in_total_idx_1, cat_active_in_xbind_idx_1 )

define BCN2¶

bn2 and cn

In [11]:

# binding & catalysis network for 2nd reaction network
# E + 2S = C

k_bind_2 = np.array([1])  # kbind = k- / k+
k_cat_2 = np.array([1])  # kcat_s2p

# bcn

total_const_idx_2 = np.array([0,])  # qe is first in totals (qe, qs, qp)
xcat_in_total_idx_2 = np.array([1, 2])  # qs is 2nd element in totals, tp is not contained.
cat_active_in_xbind_idx_2 = np.array([3])  # c_es

# binding & catalysis network for 2nd reaction network # E + 2S = C k_bind_2 = np.array([1]) # kbind = k- / k+ k_cat_2 = np.array([1]) # kcat_s2p # bcn total_const_idx_2 = np.array([0,]) # qe is first in totals (qe, qs, qp) xcat_in_total_idx_2 = np.array([1, 2]) # qs is 2nd element in totals, tp is not contained. cat_active_in_xbind_idx_2 = np.array([3]) # c_es

In [12]:

bcn2 = binding_and_catalysis(
    bn2,
    cn,
    k_bind_2,
    k_cat_2,
    np.array([total_E,]),
    total_const_idx_2,
    xcat_in_total_idx_2,
    cat_active_in_xbind_idx_2
)

bcn2 = binding_and_catalysis( bn2, cn, k_bind_2, k_cat_2, np.array([total_E,]), total_const_idx_2, xcat_in_total_idx_2, cat_active_in_xbind_idx_2 )

define BCN3¶

bn3 and cn

In [13]:

# binding & catalysis network for 2nd reaction network
# E + S = C1, S + C1 = C2

k_bind_3 = np.array([1, 1])  # kbind = k- / k+
k_cat_3 = np.array([0.1, 0.1])  # kcat_s2p
# k_cat_3 = np.array([1, 1])
# k_cat_3 = np.array([0.5, 0.5])

# bcn

total_const_idx_3 = np.array([0,])  # qe is first in totals (qe, qs, qp)
xcat_in_total_idx_3 = np.array([1, 2])  # qs is 2nd element in totals, tp is not contained.
cat_active_in_xbind_idx_3 = np.array([3, 4])  # c_es

# binding & catalysis network for 2nd reaction network # E + S = C1, S + C1 = C2 k_bind_3 = np.array([1, 1]) # kbind = k- / k+ k_cat_3 = np.array([0.1, 0.1]) # kcat_s2p # k_cat_3 = np.array([1, 1]) # k_cat_3 = np.array([0.5, 0.5]) # bcn total_const_idx_3 = np.array([0,]) # qe is first in totals (qe, qs, qp) xcat_in_total_idx_3 = np.array([1, 2]) # qs is 2nd element in totals, tp is not contained. cat_active_in_xbind_idx_3 = np.array([3, 4]) # c_es

In [14]:

bcn3 = binding_and_catalysis(
    bn3,
    cn,
    k_bind_3,
    k_cat_3,
    np.array([total_E,]),
    total_const_idx_3,
    xcat_in_total_idx_3,
    cat_active_in_xbind_idx_3
)
    

bcn3 = binding_and_catalysis( bn3, cn, k_bind_3, k_cat_3, np.array([total_E,]), total_const_idx_3, xcat_in_total_idx_3, cat_active_in_xbind_idx_3 )

Simulating the dynamic process¶

In [15]:

def get_trajs(bn, bnc, k_bind, total_S, total_P, start_time=0, end_time=20, timestamps=20):
    logxcat_init = np.log10(np.array([total_S - total_P, total_P]))  # S, P
    
    a_mat = bn.l_mat
    bnc.kbind = k_bind
    bnc.logkbind = np.log10(k_bind)
    trajs = {}
    
    time, logxcat_traj, logxbind_traj, logder_traj = bnc.get_traj(logxcat_init, start_time, end_time, timestamps, a_mat=a_mat)
    
    trajs['x_cat'] = 10 ** logxcat_traj
    trajs['x_bind'] = 10 ** logxbind_traj
    trajs['logder'] = logder_traj
    trajs['time'] = time
    
    return trajs

def get_trajs(bn, bnc, k_bind, total_S, total_P, start_time=0, end_time=20, timestamps=20): logxcat_init = np.log10(np.array([total_S - total_P, total_P])) # S, P a_mat = bn.l_mat bnc.kbind = k_bind bnc.logkbind = np.log10(k_bind) trajs = {} time, logxcat_traj, logxbind_traj, logder_traj = bnc.get_traj(logxcat_init, start_time, end_time, timestamps, a_mat=a_mat) trajs['x_cat'] = 10 ** logxcat_traj trajs['x_bind'] = 10 ** logxbind_traj trajs['logder'] = logder_traj trajs['time'] = time return trajs

In [16]:

idx_P_cn = 1

trajs_1 = get_trajs(bn1, bcn1, k_bind_1, total_S, total_P)
trajs_2 = get_trajs(bn2, bcn2, k_bind_2, total_S, total_P)
trajs_3 = get_trajs(bn3, bcn3, k_bind_3, total_S, total_P)

fig = go.Figure()

fig.add_trace(go.Scatter(x=trajs_1['time'], y=trajs_1['x_cat'][:, idx_P_cn].tolist(),
                         mode='lines', name='bn1',
                         ))

fig.add_trace(go.Scatter(x=trajs_2['time'], y=trajs_2['x_cat'][:, idx_P_cn].tolist(),
                         mode='lines', name='bn2',
                         ))

fig.add_trace(go.Scatter(x=trajs_3['time'], y=trajs_3['x_cat'][:, idx_P_cn].tolist(),
                         mode='lines', name='bn3',
                        ))

fig.update_layout(
    font=dict(
        size=27
    ),
    template='simple_white',
    width=800, height=400,
)

fig.show(renderer='notebook')

idx_P_cn = 1 trajs_1 = get_trajs(bn1, bcn1, k_bind_1, total_S, total_P) trajs_2 = get_trajs(bn2, bcn2, k_bind_2, total_S, total_P) trajs_3 = get_trajs(bn3, bcn3, k_bind_3, total_S, total_P) fig = go.Figure() fig.add_trace(go.Scatter(x=trajs_1['time'], y=trajs_1['x_cat'][:, idx_P_cn].tolist(), mode='lines', name='bn1', )) fig.add_trace(go.Scatter(x=trajs_2['time'], y=trajs_2['x_cat'][:, idx_P_cn].tolist(), mode='lines', name='bn2', )) fig.add_trace(go.Scatter(x=trajs_3['time'], y=trajs_3['x_cat'][:, idx_P_cn].tolist(), mode='lines', name='bn3', )) fig.update_layout( font=dict( size=27 ), template='simple_white', width=800, height=400, ) fig.show(renderer='notebook')

In [17]:

trajs_1['x_cat'][:, idx_P_cn].shape

trajs_1['x_cat'][:, idx_P_cn].shape

Out[17]:

(20,)

Visulize three ROP space¶

In [18]:

logmin = -6
logmax = 6
chart='x'

def sample(bn, b_vec, kcat_list=None, logmin=-6, logmax=6, npoints=100, chart='x', logvar=None):
    # randomly sample values for x
    logvar = np.random.uniform(logmin, logmax, (npoints, bn.dim_n)) if logvar is None else logvar
    # calculate log derivatives
    logders, logx = bn.logder_num(logvar, chart=chart)
    
    if isinstance(b_vec, list):
        dots_in_rop = [bn.logder_activity_num(b, logx, logders) * k for b, k in zip(b_vec, kcat_list)]
        dots_in_rop = sum(dots_in_rop)
    else:
        dots_in_rop = bn.logder_activity_num(b_vec, logx, logders)
    
    return dots_in_rop, 10 ** logx

logmin = -6 logmax = 6 chart='x' def sample(bn, b_vec, kcat_list=None, logmin=-6, logmax=6, npoints=100, chart='x', logvar=None): # randomly sample values for x logvar = np.random.uniform(logmin, logmax, (npoints, bn.dim_n)) if logvar is None else logvar # calculate log derivatives logders, logx = bn.logder_num(logvar, chart=chart) if isinstance(b_vec, list): dots_in_rop = [bn.logder_activity_num(b, logx, logders) * k for b, k in zip(b_vec, kcat_list)] dots_in_rop = sum(dots_in_rop) else: dots_in_rop = bn.logder_activity_num(b_vec, logx, logders) return dots_in_rop, 10 ** logx

In [19]:

bcn1_ROP, bcn1_x_points = sample(bn1, b_vec=np.array([0, 0, 0, 1]))
bcn2_ROP, bcn2_x_points = sample(bn2, b_vec=np.array([0, 0, 0, 1]))
bcn3_ROP, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, k_cat_3[0], k_cat_3[1]]))
bcn3_ROP_1, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, 1, 0]))
bcn3_ROP_2, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, 0, 1]))

bcn1_ROP, bcn1_x_points = sample(bn1, b_vec=np.array([0, 0, 0, 1])) bcn2_ROP, bcn2_x_points = sample(bn2, b_vec=np.array([0, 0, 0, 1])) bcn3_ROP, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, k_cat_3[0], k_cat_3[1]])) bcn3_ROP_1, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, 1, 0])) bcn3_ROP_2, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, 0, 1]))

In [20]:

E_index = 0
S_index = 1
P_index = 2

E_index = 0 S_index = 1 P_index = 2

In [21]:

ROP = go.Figure()

ROP.add_trace(go.Scatter(x=bcn1_ROP[:, S_index], y=bcn1_ROP[:, E_index],
                         mode='markers', name=r'bn1',
                         ))

ROP.add_trace(go.Scatter(x=bcn2_ROP[:, S_index], y=bcn2_ROP[:, E_index],
                         mode='markers', name=r'bn2',
                         ))

ROP.add_trace(go.Scatter(x=bcn3_ROP[:, S_index], y=bcn3_ROP[:, E_index],
                         mode='markers', name=r'bn3 combination',
                         ))

ROP.add_trace(go.Scatter(x=bcn3_ROP_1[:, S_index], y=bcn3_ROP_1[:, E_index],
                         mode='markers', name=r'bn3 C1',
                         ))

ROP.add_trace(go.Scatter(x=bcn3_ROP_2[:, S_index], y=bcn3_ROP_2[:, E_index],
                         mode='markers', name=r'bn3 C2',
                         ))


ROP.update_layout(
    font=dict(
        size=18
    ),
    template='simple_white',
    width=560, height=420,
    xaxis_title=r'$\Large{\frac{\partial \log x_C}{\partial \log q_S}}$',
    yaxis_title=r'$\Large{\frac{\partial \log x_C}{\partial \log q_E}}$',
    xaxis=dict(range=[-1, 3]),
    yaxis=dict(range=[-1, 3], scaleanchor = "x", scaleratio = 1),
)

ROP.show(renderer='notebook')

ROP = go.Figure() ROP.add_trace(go.Scatter(x=bcn1_ROP[:, S_index], y=bcn1_ROP[:, E_index], mode='markers', name=r'bn1', )) ROP.add_trace(go.Scatter(x=bcn2_ROP[:, S_index], y=bcn2_ROP[:, E_index], mode='markers', name=r'bn2', )) ROP.add_trace(go.Scatter(x=bcn3_ROP[:, S_index], y=bcn3_ROP[:, E_index], mode='markers', name=r'bn3 combination', )) ROP.add_trace(go.Scatter(x=bcn3_ROP_1[:, S_index], y=bcn3_ROP_1[:, E_index], mode='markers', name=r'bn3 C1', )) ROP.add_trace(go.Scatter(x=bcn3_ROP_2[:, S_index], y=bcn3_ROP_2[:, E_index], mode='markers', name=r'bn3 C2', )) ROP.update_layout( font=dict( size=18 ), template='simple_white', width=560, height=420, xaxis_title=r'$\Large{\frac{\partial \log x_C}{\partial \log q_S}}$', yaxis_title=r'$\Large{\frac{\partial \log x_C}{\partial \log q_E}}$', xaxis=dict(range=[-1, 3]), yaxis=dict(range=[-1, 3], scaleanchor = "x", scaleratio = 1), ) ROP.show(renderer='notebook')